> 

(N 

in 



X 



Adaptive multiresolution analysis 
based on anisotropic triangulations 

Albert Cohen, Nira Dyn, Frederic Hecht 
and Jean-Marie Mirebeau * 

January 10, 2011 



O 

(N 

C . 

, Abstract 



A simple greedy refinement procedure for the generation of data-adapted triangulations is pro- 
posed and studied. Given a function / of two variables, the algorithm produces a hierarchy of 
triangulations (I'j)j>o and piecewise polynomial approximations of / on these triangulations. The 
■ refinement procedure consists in bisecting a triangle T in a direction which is chosen so as to minimize 

the local approximation error in some prescribed norm between / and its piecewise polynomial ap- 
proximation after T is bisected. The hierarchical structure allows us to derive various approximation 
^ ■ tools such as multiresolution analysis, wavelet bases, adaptive triangulations based either on greedy 

' or optimal CART trees, as well as a simple encoding of the corresponding triangulations. We give a 

general proof of convergence in the norm of all these approximations. Numerical tests performed 
in the case of piecewise linear approximation of functions with analytic expressions or of numerical 
images illustrate the fact that the refinement procedure generates triangles with an optimal aspect 
ratio (which is dictated by the local Hessian of / in case of functions). 



1 Introduction 



Approximation by piecewise polynomial functions is a standard procedure which occurs in various ap- 
plications. In some of them such as terrain data simplification or image compression, the function to be 
approximated might be fully known, while it might be only partially known or fully unknown in other 
applications such as denoising, statistical learning or in the finite element discretization of PDE's. 

In all these applications, one usually makes the distinction between uniform and adaptive approxima- 
tion. In the uniform case, the domain of interest is decomposed into a partition where all elements have 
comparable shape and size, while these attributes are allowed to vary strongly in the adaptive case. In 
5J] ' the context of adaptive triangulations, another important distinction is between isotropic and anisotropic 

5^ \ triangulations. In the first case the triangles satisfy a condition which guarantees that they do not differ 

too much from equilateral triangles. This can either be stated in terms of a minimal value > for 
every angle, or by a uniform bound on the aspect ratio 

hx 

Pt ■= — 

rr 

of each triangle T where and r^ respectively denote the diameter of T and of its largest inscribed disc. 
In the second case, which is in the scope of the present paper, the aspect ratio is allowed to be arbitrarily 
large, i.e. long and thin triangles are allowed. In summary, adaptive and anisotropic triangulations mean 
that we do not fix any constraint on the size and shape of the triangles. 

Given a function / and a norm \\ ■ \\x of interest, we can formulate the problem of finding the optimal 
triangulation for / in the X-norm in two related forms: 
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• For a given N find a triangulation Tat with N triangles and a piecewise polynomial function (of 
some fixed degree) on Tn such that ||/ — /nWx is minimized. 

• For a given e > find a triangulation Tzv with minimal number of triangles N and a piecewise 
polynomial function fj\j such that ||/ — /atHx < £■ 

In this paper X will be the norm for some arbitrary 1 < p < oo. The exact solution to such problems is 
usually out of reach both analytically and algorithmically: even when restricting the search of the vertices 
to a finite grid, the number of possible triangulations has combinatorial complexity and an exhaustive 
search is therefore prohibitive. 

Concrete mesh generation algorithms have been developed in order to generate in reasonable time 
triangulations which are "close" to the above described optimal trade-off between error and complexity. 
They are typically governed by two intuitively desirable features: 

1. The triangulation should equidistribute the local approximation error between each triangle. This 
rationale is typically used in local mesh refinement algorithms for numerical PDE's [23]: a triangle 
is refined when the local approximation error (estimated by an a-posteriori error indicator) is large. 

2. In the case of anisotropic meshes, the local aspect ratio should in addition be optimally adapted 
to the approximated function /. In the case of piecewise linear approximation, this is achieved by 
imposing that the triangles are isotropic with respect to a distorted metric induced by the Hessian 
(Pf. We refer in particular to [7 where this task is executed using Delaunay mesh generation 
techniques. 

While these last algorithms fastly produce anisotropic meshes which are naturally adapted to the 
approximated function, they suffer from two intrinsic limitations: 

1. They are based on the evaluation of the Hessian (P f, and therefore do not in principle apply to 
arbitrary functions / £ L^'(O) for 1 < P < oo or to noisy data. 

2. They are non-hierarchical: for N > A/, the triangulation T/v is not a refinement of 71\/. 

One way to circumvent the first limitation is to regularize the function /, either by projection onto 
a finite element space or by convolution by a mollifier. However this raises the additional problem of 
appropriately tuning the amount of smoothing, in particular depending on the noise level in /. 

The need for hierarchical triangulations is critical in the construction of wavelet bases, which play an 
important role in applications to image and terrain data processing, in particular data compression 
In such applications, the multilevel structure is also of key use for the fast encoding of the information. 
Hierarchy is also useful in the design of optimally converging adaptive methods for PDE's [iT l [2T | [5l [22 ] . 
However, all these developments are so far mostly restricted to isotropic refinement methods. Let us 
mention that hierarchical and anisotropic triangulations have been investigated in [18] . yet in this work 
the triangulations are fixed in advance and therefore generally not adapted to the approximated function. 

A natural objective is therefore to design adaptive algorithmic techniques that combine hierarchy and 
anisotropy, and that apply to any function f € L^lfl), without any need for regularization. 

In this paper we propose and study a simple greedy refinement procedure that achieves this goal: start- 
ing from an initial triangulation 7)^, the procedure bisects every triangle from one of its vertices to the 
mid-point of the opposite segment. The choice of the vertex is typically the one which minimizes the new 
approximation error after bisection among the three options. 

Surprisingly, it turns out that - in the case of piecewise linear approximation - this elementary strategy 
tends to generate anisotropic triangles with optimal aspect ratio. This fact is rigorously proved in |12| 
which establishes optimal error estimates for the approximation of smooth and convex functions f G C'^, 
by adaptive triangulations T/v with N triangles. These triangulations are obtained by consecutively 
applying the refinement procedure to the triangle of maximal error. The estimates in [12] are of the form 



\\f-fN\\L.<CN-^y/\det{d^f)\\\ 
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and were already established in [TO', '4' for functions which are not necessarily convex, however based on 
triangulations which are non- hierarchical and based on the evaluation of (Pf. Note that improves 
on the estimate 

\\f-fN\\L.<CN-'\\d\f\\L., - = - + 1. (1.2) 

r p 

which can be established (see §2 in [T^) for adaptive triangulations with isotropic triangles, and which 
itselfs improves on the classical estimate 

\\f-fN\\LP<CN-'\\d^f\\L., (1.3) 

which is known to hold for uniform triangulations. 

The main objective of the present paper is to introduce the refinement procedure as well as several 
approximation methods based on it, and to study their convergence for an arbitrary function f € U' . 
In §2, we introduce notation that serves to describe the refinement procedure and define the anisotropic 
hierarchy of triangulations {'Dj)j>Q. We show how this general framework can be used to derive adaptive 
approximations of / either by triangulations based on greedy or optimal trees, or by wavelet thresholding. 
In §3, we show that as defined, the approximations produced by the refinement procedure may fail to 
converge for certain f £ and show how to modify the procedure so that convergence holds for any 
arbitrary f £ L^. We finally present in §4 some numerical tests which illustrate the optimal mesh 
adaptation, in the case of piecewise linear elements, when the refinement procedure is applied either to 
synthetic functions or to numerical images. 

2 An adaptive and anisotropic mult iresolut ion framework 
2.1 The refinement procedure 

Our refinement procedure is based on a local approximation operator At acting from LP[T) onto n,„ - the 
space of polynomials of total degree less or equal to m. Here, the parameters m > and 1 < p < oo are 
arbitrary but fixed. For a generic triangle T, we denote by (a, &, c) its edge vectors oriented in clockwise 
or anticlockwise direction so that 

a + b + c = Q. 

We define the local approximation error 

eriDp ■■= \\f - ATf\\Lv(T)- 
The most natural choice for At is the operator Bt of best U'[T) approximation which is defined by 

11/ - St/IIl^ct) = mill 11/ - 7r||ip(T). 

However this operator is non-linear and not easy to compute when p ^ 2. In practice, one prefers to use 
an operator which is easier to compute, yet nearly optimal in the sense that 

||/-^t/||lp(t) <C inf \\f-A\LHT)- (2.4) 

with C a Lebesgue constant independent of / and T . Two particularly simple admissible choices of 
approximation operators are the following: 

1. At — Pt, the (T)-orthogonal projection onto H^ , defined by Pt/ & II,„ such that Jrp{f ~ PtJ)'^ — 
for all TT g H„i. This operator has finite Lebesgue constant for all p, with C = 1 when p ~ 2 and 
C > 1 otherwise. 

2. At — It, the local interpolation operator which is defined by It/ G Hm such that /t/(7) = f{j) 
for all 7 e S {J2 ; ki e IN, = ™} where {vi,V2,V3} are the vertices of T (in the 
case m — we can take for E the barycenter of T). This operator is only defined on continuous 
functions and has Lebesgue constant C > 1 in the L°° norm. 
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All our results are simultaneously valid when At is either Pt or It (in the case where p = oo), or any 
linear operator that fulfills the continuity assumption p.4p . 

Given a target function, our refinement procedure defines by induction a hierarchy of nested trian- 
gulations {'Dj)j>o with = 2-'#('Do)- The procedure starts from the coarse triangulation Vq of f2, 

which is fixed independently of /. When O = [0, 1]^ we may split it into two symmetric triangles so that 
#(2?o) — 2. For every T £ Vj, we split T into two sub-triangles of equal area by bisection from one of 
its three vertices towards the mid-point of the opposite edge e € {a,b,c}. We denote by and the 
two resulting triangles. The choice of e e {a, 6, c} is made according to a refinement rule that selects 
this edge depending on the properties of /. We denote by TZ this refinement rule, which can therefore be 
viewed as a mapping 

n:{f,T)^e. 

We thus obtain two children of T corresponding to the choice e. is the triangulation consisting of 

all such pairs corresponding to all T £ Vj. 

In this paper, we consider refinement rules where the selected edge e minimizes a decision function 
e dT^e, f) among {a, 6, c}. We refer to such rules as greedy refinement rules. A more elaborate type of 
refinement rule is also considered in §3.3. 

The role of the decision function is to drive the generation of anisotropic triangles according to the 
local properties of /, in contrast to simpler procedures such as newest vertex bisection (i.e. split T from 
the most recently created vertex) which is independent of / and generates triangulations with isotropic 
shape constraint. 

Therefore, the choice of dT^e, f) is critical in order to obtain triangles with an optimal aspect ratio. 
The most natural choice corresponds to the optimal split 

dT(e,/) = eTi(/)^ + eT2(/)^, 

i.e. choose the edge that minimizes the resulting error after bisection. It is proved in [12] in the case 
of piecewise linear approximation, that when / is a function which is strictly convex or concave the 
refinement rule based on the decision function 

dTieJ) ^\\f - lT^if\\mT^^) + \\f - iT^fhHT^). (2.5) 

generates triangles which tend to have have an optimal aspect ratio, locally adapted to the Hessian d^f. 
This aspect ratio is independent of the norm in which one wants to minimize the error between / and 
its piecewise affine approximation. 

Remark 2.1 // the minimizer e is not unique, we may choose it among the multiple minimizers either 
randomly or according to some prescribed ordering of the edges (for example the largest coordinate pair 
of the opposite vertex in lexicographical order). 

Remark 2.2 The triangulations Vj which are generated by the greedy procedure are in general non- 
conforming, i.e. exhibit hanging nodes. This is not problematic in the present setting since we consider 
approximation in the norm which does not require global continuity of the piecewise polynomial func- 
tions. 

The refinement rule TZ defines a multiresolution framework. For a given / £ L'^iTl) and any triangle 
T we denote by 

C(T) :-{ri,T2}, 

the children of T which are the two triangles obtained by splitting T based on the prescribed decision 
function dT{e, /). We also say that T is the parent of Ti and T2 and write 

T = r{Ti) = V{T2). 

Note that 

V, :=UTeD,_,C(r). 

We also define 

V :=Uj>oPj, 
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which has the structure of an infinite binary tree. Note that Vj depends on / (except for j = 0) and on 
the refinement rule TZ, and thus V also depends on / and TZ: 

Vj^Vjif.n) and v = v{f,n). 

For notational simplicity, we sometimes omit the dependence in / and TZ when there is no possible 
ambiguity. 

2.2 Adaptive tree-based triangulations 

A first application of the multiresolution framework is the design of adaptive anisotropic triangulations 
T/v for piecewise polynomial approximation, by a greedy tree algorithm. For any finite sub-tree S G V, 
we denote by 

C{S) ■.= {TeS s.t. C{T) i 5} 
its leaves which form a partition of il. We also denote by 

I{S) ■.= S\£{S), 

its inner nodes. Note that any finite partition of by elements of V is the set of leaves of a finite sub-tree. 
One easily checks that 

#(5) = 2#{C{S)) - No. 

For each N, the greedy tree algorithm defines a finite sub-tree Sn of I? which grows from Snq '■= T^o = Tnq ) 
by adding to Sn^i the two children of the triangle T^_^ which maximizes the local U' -error eT{f)p over 
all triangles in Tzv-i- 

The adaptive partition T/v associated with the greedy algorithm is defined by 

Tn C{Sn). 

Similarly to V, the triangulation Tzv depends on / and on the refinement rule TZ, but also on p and on 
the choice of the approximation operator At- We denote by /n the piecewise polynomial approximation 
to / which is defined as Arf on each T € Tn. The global approximation error is thus given by 

11/ - /jvIIlp = ll(eT(/)p)|Up(r«)- 
Stopping criterions for the algorithm can be defined in various ways: 

• Number of triangles: stop once a prescribed N is attained. 

• Local error: stop once eT{f)p < £ for all T E Tn, for some prescribed e > 0. 

• Global error: stop once ||/ — /nWlp < £ for some prescribed e > 0. 

Remark 2.3 The role of the triangle selection based on the largest eT{f)p is to equidistrihute the local 
LP error, a feature which is desirable when we want to approximate f in LP{n) with the smallest number 
of triangles. However, it should be well understood that the refinement rule may still be chosen based 
on a decision function defined by approximation errors in norms that differ from . In particular, as 
explained earlier, the decision function (j2.5p generates triangles which tend to have have an optimal aspect 
ratio, locally adapted to the Hessian d'^f when f is strictly convex or concave, and this aspect ratio is 
independent of the norm in which one luants to minimize the error between f and its piecewise affine 
approximation. 

The greedy algorithm is one particular way of deriving an adaptive triangulation for / within the 
multiresolution framework defined by the infinite tree T). An interesting alternative is to build adaptive 
triangulations within T> which offer an optimal trade-off between error and complexity. This can be done 
when 1 < p < oo by solving the minimization problem 

mm{ eT{f)l + mS)] (2.6) 

Te£(5) 
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among all finite trees, for some fixed A > 0. In this approach, we do not directly control the number of 
triangles which depends on the penalty parameter A. However, it is immediate to see that if TV = N{X) is 
the cardinality of TJ^ = C(S*) where S* is the minimizing tree, then TJ^ minimizes the approximation 
error 

r;^:=Argmin^eT(/)^, 

where the minimum is taken among all partitions T of within V of cardinality less than or equal to A^. 

Due to the additive structure of the error term, the minimization problem (|2.6p can be performed 
in fast computational time using an optimal pruning algorithm of CART type, see [H] . In the case 
p = oo the associated minimization problem 

minj sup eT(/)oo + A#(5)|, (2.7) 

can also be solved by a similar fast algorithm. It is obvious that this method improves over the greedy 
tree algorithm: if N is the cardinality of the triangulation resulting from the minimization in (j2.6p and 
the corresponding piecewise polynomial approximation of / associated with this triangulation, we 
have 

||/-/i^||LP < Wf-fNh., 

where /jv is built by the greedy tree algorithm. 
2.3 Anisotropic wavelets 

The multiresolution framework allows us to introduce the piecewise polynomial multiresolution spaces 

Vj = Vjif,n) {g s.t. g^T e n„„ T e V^}, 
which depend on / and on the refinement rule TZ. These spaces are nested and we denote by 

y = y(/,7^) = u,>oV;■(/,7^), 

their union. For notational simplicity, we sometimes omit the dependence in / and TZ when there is no 
possible ambiguity. 

The Vj spaces may be used to construct wavelet bases, following the approach introduced in [2] and 
that we describe in our present setting. 

The space Vj is equipped with an orthonormal scaling function basis: 

^T, i^l,--- ,^(m + l)(m + 2), TeV,, 

where the for i = 1, • ■ ■ , ^{ni + l){m + 2) are supported in T and constitute an orthonormal basis 
of Urn in the sense of L^{T) for each T e P. There are several possible choices for such a basis. In the 
particular case where m = 1, a simple one is to take for T with vertices (wi,W2,W3), 

^i^{v^) = 1^1^1/2^ and (^^(v,) - -\T\-^/^VS, j + i. 
We denote by Pj the orthogonal projection onto Vj: 

TeVj i 

We next introduce for each T Vj a set of wavelets 

V-T, « = ,^("i+l)(w + 2), 

which constitutes an orthonormal basis of the orthogonal complement of Ilm{T) into Ilm{T') © Ilra{T") 
with {T',T"} the children of T. In the particular case where m = 1, a simple choice for such a basis 
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is as follows: if {vi,V2,V3) and {wi,W2,W3) denote the vertices of T' and T", with the convention that 
vi = wi and V2 = W2 denote the common vertices, the second one being the midpoint of the segment 
(^3,^3) (i.e. T has vertices (w3, W3, wi)), then 



T — 72 ' 

1 2 1,2 

,2 ._ Vj„-Vt'-Vt"+'^'t" 
Vt ■— 2 ' 

,3 yT'-yT'+yT"-yr" 

Vt ■— 2 ■ 

where cp^rp, and are the above defined scaling functions. 
The family 

^T, « = - ,^(m + l)(m + 2), TePj 

constitutes an orthonormal basis of Wj, the L^-orthogonal complement of in V^'^^. A multiscale 
orthonormal basis of Vj is given by 

{VtItsDo U {V'T}i=l,--- ,i(m+l)(m+2),- 
T£Vj,j=Q,--- ,J-1 



Letting J go to +00 we thus obtain that 



{Vt}t(^Vo U {^T}i=l,--- ,i(m+l)(m+2), 

is an orthonormal basis of the space 

v{f,n)2 :^W7Rf^''^ ^^,>^v,u,nf^''\ 

For the sake of notational simplicity, we rewrite this basis as 

(V'A)AeA, 

Note that V{f,TZ) is not necessarily dense in L'^{fl) and so V{f,Tl)2 is not always equal to L^(r2). 
Therefore, the expansion of an arbitrary function g e L^(f2) in the above wavelet basis does not always 
converge towards g in L^(r2). The same remark holds for the convergence of the wavelet expansion of 
an arbitrary function g E L'P{n) (or C{n) in the case p = 00): L^-convergence holds when the space 



coincides with LP{Q) (or contains C(il) in the case p = 00), since we have 

/-V dx^Px =\\f-P,f\\L.<C M ||/-g||Lp, 

— ' LP geVj 



with C the Lebesgue constant in (|2.4p for the orthogonal projector. 

A sufficient condition for such a property to hold is obviously that the size of all triangles goes to 
as the level j increases, i.e. 

lim sup diam(T) = 0. 

j->-+oo TeVj 

However, this condition might not hold for the hierarchy (2?j)j>o produced by the refinement procedure. 
On the other hand, the multiresolution approximation being intrinsically adapted to /, a more reasonable 
requirement is that the expansion of / converges towards / in LP{rt) when / e Lp{^) (or C(r2) in the 
case p = 00). This is equivalent to the property 



/ey(/,7^) 



p. 
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We may then define an adaptive approximation of / by thresliolding its coefficients at some level e > 0: 

l/Al>e 

where fx '■= (/jV'a)- When measuring the error in the norm, a more natural choice is to perform 
thresholding on the component of the expansion measured in this norm, defining therefore 

fe ■■= Y Mx- 

\\Mx\\lp>s 

Wc shall next see that the condition / € V{f, TZ)p also ensures the convergence of the tree-based adaptive 
approximations /jv and towards / in Lp{Q). We shall also see that this condition may not hold for 
certain functions /, but that this difficulty can be circumvented by a modification of the refinement 
procedure. 



3 Convergence analysis 
3.1 A convergence criterion 

The following result relates the convergence towards / of its approximations by projection onto the 
spaces Vj{f,TZ), greedy and optimal tree algorithms, and wavelet thresholding. This result is valid for 
any refinement rule TZ. 

Theorem 3.1 Let TZ : (/, T) i— > e be an arbitrary refinement rule and let f G L'''{il). The following 
statements are equivalent: 

(i) /e^/(/,7^)p. 

(ii) The greedy tree approxim,ation converges: limjv->.+oo 11/ ~ fwllLP = 0. 

(iii) The optimal tree approximation converges: limjv->-i-cio ||/ ~ /jvlli''' ~ 0- 
In the case p = 2, they are also equivalent to: 

(iv) The thresholding approxim,ation converges: lime^.o ||/~ /elU^ =0. 

Proof: Clearly, (ii) implies (iii) since ||/ — /jvlUp ^ 11/ ~ fNWhp- Since the triangulation 2?jv is a 
refinement of 7^, we also find that (ni) implies (i) as infggvw 11/ ~ slli" ^ 11/ ^ //vlU''- 
We next show that (i) implies (ii). We first note that a consequence of (i) is that 

lim sup eT(/)p = 0. 

It follows that for any Jj > 0, there exists N{ri) such that for A'' > N{r]), all triangles T e Tn satisfy 

eT(/)p < V- 

On the other hand, (i) means that for all e > 0, there exists J = J{e) such that 

inf \\f-g\\L.<e. 

For A'" > N{ri), we now split Tn into TJ^ U TJ^ where 

r+ i^Tn^ (Ujyj'Dj) and Tj^ := Tn n {Uj<jVj). 
We then estimate the error of the greedy algorithm by 

\\f-fN\\i. = E eT(/)^+ E 

Ter+ tgTj^ 
Ter+ 

< CP inf \\f-g\\l^+rjPi^{r^) 
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where C is the stabihty constant of <\2A\ . This imphes (ii) since for any (5 > 0, we can first choose e > 
such that C^eP < S/2, and then choose 77 > such that 2'^ '■''^ NqtjP < S/2. When p = 00 the estimate is 
modified into 

\\f ~ fwh^ <max{Ce,77}, 

which also implies (ii) by a similar reasoning. 

We finally prove the equivalence between (i) and (iv) when p — 2. Property (i) is equivalent to the 
convergence of the orthogonal projection Pjf to f as j +00, or equivalently of the partial sum 

\M<j 

where | A| stands for the scale level of the wavelet V'a- Since ('0a)a6A is an orthonormal basis of V{f, 71)2, 
the summability and limit of X^agA /a^a do not depend on the order of the terms. Therefore (i) is 
equivalent to the convergence of to /. o 

Remark 3.2 The equivalence between statements (i) and (iv) can be extended to 1 < p < 00 by showing 
that {')jjx)xeA is an -unconditional system. Recall that this property means that there exists an absolute 
constant C > such that for any finitely supported sequences (ca) and (dx) such \cx\ < \dx\ for all X, 
one has 

II ^ca^AaIIlp < C\\ y^^dxipx\\Lp- 
A consequence of this property is that if f Cz can be expressed as the limit 

/ = lim V fxipx, 

7^4-00 ^ ^ 

|A|<J 

then any rearrangement of the series X) /aV'a converges towards f in LP. This easily implies the equiv- 
alence between (i) and (iv). The fact that (V'a)agA is an unconditional system is well known for Haar 
systems U9f which correspond to the case m — 0, and can be extended to m > in a straightforward 
manner. 



3.2 A case of non-convergence 

We now show that if we use a greedy refinement rule TZ based on a decision function either based 
on the interpolation or projection error after bisection, there exists functions / G LP(fl) such that 
/ ^ V{f,TZ)p. Without loss of generality, it is enough to construct / on the reference triangle Trof of 
vertices {(0,0), (1,0), (1, 1)}, since our construction can be adapted to any triangle by an afhne change 
of variables. 

Consider first a decision function defined from the interpolation error after bisection, such as (|2.5p . 
or more generally 

dTif,e) 11/ - /ti/IIL(t,i) + 11/ ^ Ir^fWUw- 

Let / be a continuous function which is not identically on Tj-ef and which vanishes at all points {x, y) such 
that X = 2^ for fc = 0, 1, • • • , 2m, where m is the degree of polynomial approximation. For such an /, it is 
easy to see that Ixf ~ and that /t'/ — for all subtriangles T' obtained by one bisection of Trof. This 
shows that there is no preferred bisection. Assuming that we bisect from the vertex (0, 0) to the opposite 
mid-point (1, 5), we find that a similar situation occurs when splitting the two subtriangles. Iterating 
this observation, we see that an admissible choice of bisections leads after j steps to a triangulation T>j of 
Tref consisting of the triangles Tj^k with vertices {(0, 0), (1, 2-^k), (1, 2-^{k + 1))} with fc = 0, • • • , 2-' - 1. 
On each of these triangles / is interpolated by the null function and therefore by (|2.4p the best L°° 
approximation in Vj does not converge to / as j — > +00, i.e. / ^ V{f, TZ)oo- It can also easily be checked 

that y(/,7^)p. 

Similar counter-examples can be constructed when the decision function is defined from the pro- 
jection error, and has the form 

rfT(/,e) ||/-Pti/|IL(ti) + II/-^t,^/IIL(t,=)- 
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Here, wc describe such a construction in the case m = 1. We define / on 7?. as a function of the first 
variable given by 



f{x,y) = u{x), if X G 



, f{x,y) 



if X e 



2' 



where u is a non-trivial function in L^([0, ^]) such that 



u{x)dx 



xu(x)dx 



x^u{x)dx = 0. 



A possible choice is u{x) = L^iAx — 1) = 160x^ — 120x'^ + 24.T — 1 where L3 is the Legendre polynomial 
of degree 3 defined on [—1,1]. With this choice, we have the following result. 

Lemma 3.3 Let T be any triangle such that its vertices have x coordinates either (0, 5,1) or (0, 1, 1) or 
(i, 1, 1). Then f is orthogonal to Hi in L'^{T). 



Proof: Define 



Po := |(x, 



y)eT, xG 



and Ti 



:=|(x, 



y)eT, xG ( -,1 



Then, with v{x,y) being either the function 1 or x or y, we have 

It /(^' y)v{x, y)dxdy = Jj,^ u{x)v{x, y)dxdy + /^^ u{x)v{x, y)dxdy 



= u{x)qo{x)dx + u{x)qi{x)dx, 



with 



90 



(x) := / V 



{x,y)dy, qi{x) :-- 



v{x,y)dy, 



where Ti^^ = {y ■ (x, y) G T,} for i = 0,1. The functions qo and qi are polynomials of degree at most 2 
and we thus obtain from the properties of u that Jj,fv = 0. o 



The above lemma shows that for any of the three possible choices of bisection of Trcf based on the 
decision function, the error is left unchanged since the projection of / on all possible sub-triangle is 0. 
There is therefore no preferred choice, and assuming that we bisect from the vertex (0, 0) to the opposite 
mid-point (1,5), then we see that a similar situation occurs when splitting the two subtriangles. The 
rest of the arguments showing that / ^ V{f, TZ)p are the same as in the previous counter-example. 

The above two examples of non-convergence reflect the fact that when / has some oscillations, the 
refinement procedure cannot determine the most appropriate bisection. In order to circumvent this 
difficulty one needs to modify the refinement rule. 



3.3 A modified refinement rule 

Our modification consists of bisecting from the most recently generated vertex of T, in case the local 
error is not reduced enough by all three bisections. More precisely, we modify the choice of the bisection 
of any T as follows: 

Let e be the edge which minimizes the decision function drie, /). If 

{eT.{f)l + er^{fXy^'<9eTif), 

we bisect T towards the edge e (greedy bisection). Otherwise, we bisect T from its most recently generated 
vertex (newest vertex bisection). Here ^? is a fixed number in (0, 1). In the casep = 00 we use the condition 

max{eTi(/),eT2(/)} < Oerif). 
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This new refinement rule benefits from the mesh size reduction properties of newest vertex bisection. 
Indeed, a property illustrated in Figure 1 is that a sequence {BNN} of one arbitrary bisection (B) 
followed by two newest vertex bisections (N) produces triangles with diameters bounded by half the 
diameter of the initial triangle. A more general property - which proof is elementary yet tedious - is 
the following: a sequence of the type {BNB ■ ■ ■ BN} of length k + 2 with a newest vertex bisection at 
iteration 2 and fc + 2 produces triangles with diameter bounded by (1 — 2^'^') times the diameter of the 
initial triangle, the worst case being illustrated in Figure 2 with fc = 3. 




Figure 1: diameter reduction by a {BNN} sequence 




Figure 2: diameter reduction by a {BNBBN} sequence 
(the dark triangle has diameter at most 7/8 times the initial diameter) 



Our next result shows that the modified algorithm now converges for any f £ L^. 
Theorem 3.4 With TZ defined as the modified bisection rule, we have 

feV{f,n)p, 

for any f e LP(f2) ( or C{Q) when p — oo). 

Proof: We first give the proof when p < oo. For each triangle T S T>j with j > 1, we introduce the two 
quantities 

(^(T) — . .^p . ff.p and f3{T) := ^ — ^ , 

eT(/)p + eT'(/)p e-p(T)(/)p 

where T' is the "brother" of T, i.e. C{V{T)) = {r,T'}. When a greedy bisection occurs in the split of 
■p(T), we have P{T) < 9^. When a newest vertex bisection occurs, we have 

ini^en„ ||/ - T^WLPiviT)) 



where C is the constant of 

We now consider a given level index j > and the triangulation Vj. For each T E Vj, wc consider 
the chain of nested triangles (r„)^^Q with Tj — T and r„_i — V{Tn), n = j, j — 1, • • • ,1. We define 

j j 
0a(T) = W a{Tn) and 0/3(r) = [] /3(r„). 

n—l n—1 

It is easy to see that 

eT(f)P 

0a{T)0(5{T) = 

eTo(/)p 

so that 

eriffp < Co0a(T)0/3(T), 
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with Co '■— maxj-Qg-p,-, eTo{f)^. It is also easy to check by induction on j that 

0c,{T) = 0a{T) = ■■■= Y MT) - #I?o. 

We denote by fj the approximation to / in Vj defined by fj = Arf on all T E Tj so that 

In order to prove that fj converges to / in i^, it is sufficient to show that the sequence 

Ej := mEix min{0/3(r), diam(T)}, 

tends to as j grows. Indeed, if this holds, we split Vj into two sets Vj' and VJ over which 0f3(T) < ej 
and diam(r) < Ej respectively. We can then write 

\\.f-m% = E ^^(/)p+ E ^Tifr, 

Tev+ Tevr 
< Cos, Y MT) + CP Y JlJI-^-^llinT) 



< Come, + CP Y ll/-^lli.m' 



where C is the constant of ()2.4p . Clearly the first term tends to and so does the second term by standard 
properties of spaces since the diameter of the triangles in VJ goes to 0. It thus remains to prove that 



lim Ej — 0. 



Again we consider the chain (T'„)^^q which terminates at T, and we associate to it a chain (qn)l=o where 
g„ = 1 or 2 if bisection of r„ is greedy or newest vertex respectively. If r is the total number of 2 in the 
chain (q„) we have 

0/3(r) < cp''ep(J-''\ 



with C the constant in (|2.4[) . Let a fc > be a fixed number, large enough such that 

ce''-^ < 1. 

We thus have 

0P{T) < (CP0P'-^-^^Y0P'-^'''''^ < 0P<-J-'-''K (3.8) 

We now denote by I the maximal number of disjoint sub-chains of the type {ui, 2, 1/2, t'a, • • • , I'q, 2) with 
Vj e {1, 2} and of length q + 2 < 2k + 3 which can be extracted from ((j'n)^=o- Fi'oni the remarks on the 
diameter reduction properties of newest vertex bisection, we see that 

diam(T) < B{1 - 2-'^^)\ 

with B :— maxj'Qgx'o diam(To) a fixed constant. On the other hand, it is not difficult to check that 

r <3l + 3 + ll-!l, (3.9) 
Zk 

Indeed let be the total number of 1 in the sequence ((7„) which are not preceeded by a 2, and let 
be the size of the series of 1 following the z-th occurence of 2 in (qn) for i — 1, - ■ ■ , r. Note that some 
might be 0. Clearly we have 

j = ao + ai + ■ ■ ■ + ar + r. 
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From the above equality, the number of i such that > 2fc is less than and therefore there is at least 
m > r — indices {io, • • • , im-i} such that Ui < 2k. Denoting by Pi the position of the i-th occurence 
of 2 (so that Pi+i — I3i + ai + I), we now consider the disjoint sequences of indices 

St = {/3i3t, ■ ■ ■ jft3t + 2}j t — ' ' ' 

There is at least y — 1 such sequences within {1, • • • , j} and by construction each of them contains a 
sequence of the type (j^i, 2, 1/2, t's, • ■ ■ , i^q, 2) with i/j G {1, 2} and of length g + 2 < 2A; + 3. Therefore the 
maximal number of disjoint sequences of such type satisfies 

- 3^ 2k ^ 
which is equivalent to p.9|) . Therefore according to p.8p 

0l3{T) < < 6»p(5-3(/+i)fc+i) 

If 3(Z + l)fc< f, we have 

0/3(r) < e'^. 

On the other hand, if 3{l + l)k > |, we have 

diam(T) < B{1 - 2-^'')^-\ 
We therefore conclude that Sj goes to as j grows, which proves the result for p < 00. 
We briefly sketch the proof for p = 00, which is simpler. We now define /3(T) as 



e-v{T){f)c 



so that li{T) < 6* if a greedy bisection occurs in the split of V{T). With the same definition of 0(3 {T) we 
now have 

erifU < Co0/3(r), 

where Co '■— maxT^gp^ exaif) co- With the same definition of ej and splitting of 2?j, we now reach 

11/ - < max < CoEj, C max inf \\ f - tt\\ {t) ? , 
I Tevr iren„, I 

which again tends to if ej tends to and / is continuous. The proof that Sj tends to as j grows is 
then similar to the case p < 00. o 

Remark 3.5 The choice of the parameter 9 < 1 deserves some attention: if it is chosen too small, then 
most bisections are of type N and we end up with an isotropic triangulation. In the case m = 1, a 
proper choice can he found by observing that when f has smoothness, it can he locally approximated 
by a quadratic polynomial g G 112 with eT{f)p ~ eT{<l)p when the triangle T is small enough. For such 
quadratic functions q S 112, one can explicitely study the minimal error reduction which is always ensured 
by the greedy refinement rule defined a given decision function. In the particular case p = 2 and with 
the choice At = Pt which is considered in the numerical experiments of explicit formulas for the 
error \\q — Pt9||l2(t) can he obtained by formal computing and can be used to prove a guaranteed error 
reduction by a factor 0* = |. It is therefore natural to choose 9 such that 9* < 9 < 1 (for example 9 — ^) 
which ensures that bisections of type N only occur in the early steps of the algorithm, when the function 
still exhibits too many oscillations on some triangles. 
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4 Numerical illustrations 



The following numerical experiments were conducted with piecewise linear approximation in the 
setting: we use the L^-based decision function 

dT{f,e) := \\f - PT^fWl^T^-^ + Wf - Pr^fWh^T^y 

and we take for At the L^(T)-orthogonal projection Pt onto Hi. In these experiments, the function / is 
either a quadratic polynomial or a function with a simple analytic expression which allows us to compute 
the quantities eT(/)2 and dT{e, f) without any quadrature error, or a numerical image in which case the 
computation of these quantities is discretized on the pixel grid. 

4.1 Quadratic functions 

Our first goal is to illustrate numerically the optimal adaptation properties of the refinement procedure 
in terms of triangle shape. For this purpose, we take / = q a quadratic form i.e. an homogeneous 
polynomial of degree 2. In this case, all triangles should have the same aspect ratio since the Hessian is 
constant. In order to measure the quality of the shape of a triangle T in relation to q, we introduce the 
following quantity: if (a, b, c) are the edge vectors of T, we define 

frp^ _ max{|q(a)|,|q(6)|,|q(c)|} 

''^^^•^ \T\VWm ' 

where det(q) is the determinant of the 2x2 symmetric matrix Q associated with q, i.e. such that 

q(u) = {Qu, u) 

for all u G M^. Using the reference triangle and an affine change of variables, it is proved in §2 of [T^] 
that 

eT(q)p ^ |rr+^Pq(T)V|det(q)|, 

with equivalence constants independent of q and T. Therefore, if T is a triangle of given area, its shape 
should be designed in order to minimize pq{T). 

In the case where q is positive definite or negative definite, Pq{T) takes small values when T is 
isotropic with respect to the metric \{x^y)\q := •\/|q(a;, y)\, the minimal value being attained for an 

equilateral triangle for this metric. Specifically, we choose q(a;, y) :— + lOOj/^ and display in Figure 
3 (left) the triangulation obtained after j = 8 iterations of the refinement procedure, starting with 
a triangle which is equilateral for the euclidean metric (and therefore not adapted to q). Triangles such 
that Pq{T) < 4-\/3 (at most 3 times the minimal value) are displayed in white, others in grey. We observe 
that most triangles produced by the refinement procedure are of the first type and therefore have a good 
aspect ratio. 



Figure 3: 2?8 for q{x,y) :— + lOOy (left) and q{x,y) := x — lOy (right) 



The case of a quadratic function of mixed signature is illustrated in Figure 3 (right) with q(a;,y) := 

|q| 



— lOy^. For such quadratic functions, triangles which are isotropic with respect to the metric 
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have a low value of pq, where |q| denotes the positive quadratic form associated to the absolute value \Q\ 
of the symmetric matrix Q associated to q. Recall for any symmetric matrix Q there exists Ai,A2 € K 
and a rotation R such that 

and the absolute value |Q| is defined as 

101-^" Co' |A°|)«' 

In the present case, R — I and |q|(x,y) — x'^ + lOy^. 

But one can also check that pq is left invariant by any linear transformation with eigenvalues {t, j) for 
any t > and eigenvectors {u, v) such that q(w) — q(w) = 0, i.e. belonging to the null cone of q. More 
precisely, for any such transformation ip and any triangle T, one has pq('0(T)) = pq{T). In our example 
we have u = (VTO, 1) and v = (\/T0, —I))- Therefore long and thin triangles which are aligned with these 
vectors also have a low value of pq. Triangles T such that p|q|(T) < 4-\/3 are displayed in white, those 
such that Pq(T) < 4a/3 while p|q|(T) > 4-\/3 - i.e. adapted to q but not to |q| - are displayed in grey, 
and the others in dark. We observe that all the triangles triangles produced by the refinement procedure 
except one are either of the first or second type and therefore have a good aspect ratio. 

4.2 Sharp transition 

We next study the adaptive triangulations produced by the greedy tree algorithm for a function / dis- 
playing a sharp transition along a curved edge. Specifically we take 

f{x, y) = fs{x, y) := gs{\/ x"^ +y^), 

where gs is defined by g5{r) = ^^j- for < r < 1, ^^(l + 5 -\- r) = —^niliill. fQj- r > 0, and gs is a 
polynomial of degree 5 on [1, 1 + (5] which is determined by imposing that gs is globally . The parameter 
5 therefore measures the sharpness of the transition as illustrated in Figure 4. It can be shown that the 
Hessian of fs is negative definite for \J x"^ + y"^ < 1+5/2, and of mixed type for 1+5/2 < ^ x^ + y"^ < 2+S. 



Figure 4: The function gs, for 5 = 0.02, 0.03, 0.07, 0.2. 

Figure 5 displays the triangulation Tioooo obtained after 10000 steps of the algorithm for 5 = 0.2. In 
particular, triangles T such that Pq{T) < 4 - where q is the quadratic form associated with cP f measured 
at the barycenter of T - are displayed in white, others in grey. As expected, most triangles are of the first 
type and therefore well adapted to /. We also display on this figure the adaptive isotropic triangulation 
produced by the greedy tree algorithm based on newest vertex bisection for the same number of triangles. 

Since / is a function, approximations by uniform, adaptive isotropic and adaptive anisotropic 
triangulations all yield the convergence rate 0(iV^^). However the constant 

C := limsupTVll/- /Ar||L2, 

strongly differs depending on the algorithm and on the sharpness of the transition, as illustrated in 
the table below. We denote by Cu, Ci and Ca the empirical constants (estimated by iV||/ — /Ar||2 
for N — 8192) in the uniform, adaptive isotropic and adaptive anisotropic case respectively, and by 
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Figure 5: Tioooo (left), detail (center), isotropic triangulation (right). 



U{f) 11^2/11^2, /(/) := |l(i^/|li2/3 and A{f) := || A/|det(d2/)| 11^2/3 the theoretical constants suggested 
by (|1.3I) . (II. 2p and We observe that Cu and C/ grow in a similar way as U{f) and /(/) as 5 ^ 

(a detailed computation shows that U{f) ~ 10.37 and /(/) « 14.01(5"^/^). In contrast C^i and 
A{f) remain uniformly bounded, a fact which reflects the superiority of the anisotropic mesh as the layer 
becomes thinner. 
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U{f) 


Hf) 


Mf) 


Cu 


Ci 


Ca 


0.2 


103 


27 


6.75 


7.87 


1.78 


0.74 


0.1 


602 


60 


8.50 


23.7 


2.98 


0.92 


0.05 


1705 


82 


8.48 


65.5 


4.13 


0.92 


0.02 


3670 


105 


8.47 


200 


6.60 


0.92 



4.3 Numerical images 

We finally apply the greedy tree algorithm to numerical images. In this case the data / has the form 
of a discrete array of pixels, and the (T)-orthogonal projection is replaced by the £^(S'7-)-orthogonal 
projection, where St is the set of pixels with centers contained in T. The approximated 512 x 512 image 
is displayed in Figure 6 which also shows its approximation Jn by the greedy tree algorithm based on 
newest vertex bisection with N — 2000 triangles. The systematic use of isotropic triangles results in 
strong ringing artifacts near the edges. 




Figure 6: The image "peppers" (left), /2000 with newest vertex (right). 



We display in Figure 7 the result of the same algorithm now based on our greedy bisection procedure 
with the same number of triangles. As expected, the edges are better approximated due to the presence of 
well oriented anisotropic triangles. Yet artifacts persist on certain edges due to oscillatory features in the 
image which tend to mislead the algorithm in its search for triangles with good aspect ratio, as explained 
in §3.2. These artifacts tend to disappear if we use the modified refinement rule proposed in §3.3 as also 
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Figure 7: /2000 with greedy bisection (left), modified procedure (right). 



illustrated on Figure 7. This modification is thus useful in the practical application of the algorithm, in 
addition of being necessary for proving convergence convergence of the approximations towards any 
function. Note that encoding a triangulation resulting from N iterations of the anisotropic refinement 
algorithm is more costly than for the newest vertex rule: the algorithm encounters at most 2N triangles 
and for each of them, one needs to encode one out of four options (bisect towards edge a or 6 or c or 
not bisect), therefore resulting into AN bits, while only two options need to be encoded when using the 
newest vertex rule (bisect or not), therefore resulting into 2N bits. In the perspective of applications to 
image compression, another issue is the quantization and encoding of the piecewise affine function as well 
as the treatment of the triangular visual artifacts that are inherent to the use of discontinuous piecewise 
polynomials on triangulated domains. These issues will be discussed in a further work specifically dealing 
with image applications. 

5 Conclusions and perspectives 

In this paper, we have studied a simple greedy refinement procedure which generates triangles that tend to 
have an optimal aspect ratio. This fact is rigorously proved in |12j , together with the optimal convergence 
estimate for the adaptive triangulations constructed by the greedy tree algorithm in the case where 
the approximated function / is and convex. Our numerical results illustrate these properties. 

In the present paper we also show that for a general / S the refinement procedure can be misled by 
oscillations in /, and that this drawback may be circumvented by a simple modification of the refinement 
procedure. This modification appears to be useful in image processing applications, as shown by our 
numerical results. 

Let us finally mention several perspectives that are raised from our work, and that are the object of 
current investigation: 

1. Conforming triangulations: our algorithm inherently generates hanging nodes, which might not 
be desirable in certain applications, such as numerical discretization of PDE's where anisotropic 
elements are sometimes used [3] . When using the greedy tree algorithm, an obvious way of avoiding 
this phenomenon is to bisect the chosen triangle together with an adjacent triangle in order to 
preserve conformity. However, it is no more clear that this strategy generates optimal triangulations. 
In fact, we observed that many inappropriately oriented triangles can be generated by this approach. 
An alternative strategy is to apply the non-conforming greedy tree algorithm until a prescribed 
accuracy is met, followed by an additional refinement procedure in order to remove hanging nodes. 

2. Discretization and encoding: our work is in part motivated by applications to image and terrain 
data processing and compression. In such applications the data to be approximated is usually given 
in discrete form (pixels or point clouds) and the algorithm can be adapted to such data, as shown 
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in our numerical image examples. Key issues which need to be dealt with are then (i) the efficient 
encoding of the approximations and of the triangulations using the tree structure in a similar spirit 
as in [11 and (ii) the removal of the triangular visual artifacts due to discontinuous piecewise 
polynomial approximation by an appropriate post-processing step. 

3. Adaptation to curved edges: one of the motivation for the use of anisotropic triangulations is the 
approximation of functions with jump discontinuities along an edge. For simple functions, such as 
characteristic functions of domains with smooth boundaries, the L^-error rate with an optimally 
adapted triangulation of N elements is known to be 0{N~p). This rate reflects an C(l) error 
concentrated on a strip of area 0{N~'^) separating the curved edge from a polygonal line. Our 
first investigations in this direction indicate that the greedy tree algorithm based on our refinement 
procedure cannot achieve this rate, due to the fact that bisection does not offer enough geometrical 
adaptation. This is in contrast with other splitting procedures, such as in |13| in which the direction 
of the new cutting edge is optimized within an infinite range of possible choices, or |15| where the 
number of choices grows together with the resolution level. An interesting question, addressed 
and partially answered in §9.1 of [24], is thus to understand if the optimal rate for edges can be 
achieved by a splitting procedure with a small and fixed number of choices similar to our refinement 
procedure, which would be beneficial from both a computational and encoding viewpoint. 
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